In this workbook we shall explore a process for selecting GMPEs for application to Italy following these steps:
Pre-selection of the GMPEs (no notebook required)
Comparison of the GMPEs for 'typical' scenarios
Comparison of the GMPEs against data
In [ ]:
# Setup
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
# OpenQuake tools
from openquake.hazardlib.imt import from_string
from openquake.hazardlib.gsim import get_available_gsims
# Import the tools for building a database from a flatfile
import smtk.sm_database_builder as sdb
from smtk.parsers.general_flatfile_parser import GeneralFlatfileParser
from smtk.database_visualiser import db_magnitude_distance
# Import the selection tools
from smtk.strong_motion_selector import SMRecordSelector, rank_sites_by_record_count
# Import the residual tools
import smtk.residuals.gmpe_residuals as res
import smtk.residuals.residual_plotter as rspl
# Import the trellis plotting tools
import smtk.trellis.configure as rcfg
import smtk.trellis.trellis_plots as trpl
#As a little digression, if you want to have more control over the colour/marker cycle then this
# snippit of code can help
from cycler import cycler
import matplotlib
matplotlib.rcParams["axes.prop_cycle"] = \
cycler(u'color', ['b', 'r', 'g', "c", "m", "y", "k", # Order of colours
'b', 'r', 'g', "c", "m", "y", "k"]) +\
cycler(u'linestyle', ["-", "-", "-", "-", "-", "-", "-", # Order of markers
"--", "--", "--", "--", "--", "--", "--"])
Take a look at the GMPEs in OpenQuake and thei tectonic region type
In [ ]:
available_gsims = get_available_gsims()
for gsim_name in available_gsims:
if gsim_name != "GMPETable":
gsim = available_gsims[gsim_name]()
print("%s - %s" % (gsim_name, gsim.DEFINED_FOR_TECTONIC_REGION_TYPE))
... OK that's a lot of GMPEs. What about just the active shallow crust?
In [ ]:
for gsim_name in available_gsims:
if gsim_name != "GMPETable":
gsim = available_gsims[gsim_name]()
if "Active Shallow Crust" in gsim.DEFINED_FOR_TECTONIC_REGION_TYPE:
print("%s - %s" % (gsim_name, gsim.DEFINED_FOR_TECTONIC_REGION_TYPE))
Setup the rupture
In [ ]:
rupture1 = rcfg.GSIMRupture(magnitude=6.5,
dip=60,
aspect=1.2,
rake=-90.,
ztor=0.0,
hypocentre_location=(0.5, 0.5))
# Define the target points as a line perpendicular to the strike of the fault
rupture1.get_target_sites_line(250., 1.0, vs30=800., line_azimuth=90.)
rupture1.plot_model_configuration()
Choose the GMPE List, Intensity Measure Types and Spectral Periods
In [ ]:
gsim_list = ["AkkarEtAlRjb2014",
"CauzziEtAl2014",
"BindiEtAl2014Rjb",
"BooreEtAl2014",
"ChiouYoungs2014",
"KothaEtAl2016Italy",
"AbrahamsonEtAl2014",
"CampbellBozorgnia2014"]
imt_list = ["PGA", "SA(0.1)", "SA(0.2)", "SA(0.3)", "SA(0.5)", "SA(1.0)", "SA(2.0)"]
periods = [0.05, 0.075, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19,
0.20, 0.22, 0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.36, 0.38, 0.40, 0.42, 0.44, 0.46, 0.48, 0.5,
0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
1.9, 2.0, 3.0, 4.0, 5.0, 7.5, 10.0]
In [ ]:
trpl.DistanceIMTTrellis.from_rupture_model(rupture1, gsim_list, imt_list, figure_size=(10,10))
Plot the Attenuation with Distance
In [ ]:
rupture2 = rcfg.GSIMRupture(magnitude=6.5,
dip=60,
aspect=1.2,
rake=-90.,
ztor=0.0,
hypocentre_location=(0.5, 0.5))
# Define the target points as a line perpendicular to the strike of the fault
rupture2.get_target_sites_point(5., "rrup", vs30=800., line_azimuth=90.)
rupture2.plot_model_configuration()
Plot the Scenario Spectrum
In [ ]:
rupture2 = rcfg.GSIMRupture(magnitude=6.5,
dip=60,
aspect=1.2,
rake=-90.,
ztor=0.0,
hypocentre_location=(0.5, 0.5))
# Define the target points as a line perpendicular to the strike of the fault
rupture2.get_target_sites_point(5., "rrup", vs30=800., line_azimuth=270.)
rupture2.plot_model_configuration()
trpl.MagnitudeDistanceSpectraTrellis.from_rupture_model(rupture2, gsim_list, periods, distance_type="rrup")
In [ ]:
# Setup the rupture
rupture2 = rcfg.GSIMRupture(magnitude=6.5,
dip=60,
aspect=1.2,
rake=-90.,
ztor=0.0,
hypocentre_location=(0.5, 0.5))
# Define the target points as a line perpendicular to the strike of the fault
_ = rupture2.get_target_sites_point(5., "rrup", vs30=800., line_azimuth=90.)
# Get all of the parameters the rupture needs for the calculation (i.e. site, distance, rupture)
sctx, rctx, dctx = rupture2.get_gsim_contexts()
# Merge the site information into a dictionary
sctx.__dict__.update(rctx.__dict__)
for val in dctx.__dict__:
if getattr(dctx, val) is not None:
setattr(dctx, val, getattr(dctx, val)[0])
# Magnitudes from 4.5 to 8.0 in increments of 0.1
magnitudes = np.arange(4.5, 8.1, 0.1)
trpl.MagnitudeIMTTrellis(magnitudes, dctx.__dict__,
gsim_list, imt_list, sctx.__dict__,
figure_size=(12, 12))
In [ ]:
rupture3 = rcfg.GSIMRupture(magnitude=6.0,
dip=40,
aspect=1.2,
rake=90.,
ztor=2.0,
hypocentre_location=(0.5, 0.5))
# Define the target points as a line perpendicular to the strike of the fault
rupture3.get_target_sites_line(250., 1.0, vs30=800., line_azimuth=90.)
rupture3.plot_model_configuration()
In [ ]:
trpl.DistanceIMTTrellis.from_rupture_model(rupture3, gsim_list, imt_list, figure_size=(10,10))
In [ ]:
rupture4 = rcfg.GSIMRupture(magnitude=6.0,
dip=40,
aspect=1.2,
rake=90.,
ztor=2.0,
hypocentre_location=(0.5, 0.5))
# Define the target points as a line perpendicular to the strike of the fault
rupture4.get_target_sites_point(7., "rrup", vs30=800., line_azimuth=90.)
rupture4.plot_model_configuration()
In [ ]:
trpl.MagnitudeDistanceSpectraTrellis.from_rupture_model(rupture4, gsim_list, periods, distance_type="rrup")
Look at the database file core_flatfile_ngawest2.csv.
This file is the most comprehensive format, so there are many columns that are included but are not critical.
For comparison against GMPEs some are:
In [ ]:
# Path to the file
flatfile = "./data/core_flatfile_ngawest2.csv"
# Setup the builder - needs the type of parser and then the path to store the database
builder = sdb.SMDatabaseBuilder(GeneralFlatfileParser, "./data/ngawest_db")
# Parse the metadata - needs an "ID", "A basic name" and the path where to find the flatfile
builder.build_database("NGAWEST2", "Basic NGA West 2 Flatfile ", flatfile)
In [ ]:
# Parse the spectra - needs to know the component (e.g. "Geometric", "GMRotD50", "SARotD50")
builder.build_spectra_from_flatfile("SARotD50", damping="05", units="g")
Not all the records in the database may be good for use. We can initially filter out some that may not be suitable.
Here we use the strong motion selector and apply two filterings:
Only use records with $100 \leq V_{S30} \left( {m/s} \right) \leq 1500$
Only use records with $0.1 \leq R_{RUP} \left( {km} \right) \leq 500$
In [ ]:
# Now load in the databse
import cPickle
db1 = cPickle.load(open("./data/ngawest_db/metadatafile.pkl", "r"))
print("Initial Database has %g records" % len(db1))
# Eliminate records with missing Vs30 or excessively low or high values
selector = SMRecordSelector(db1)
db1.records = selector.select_within_vs30_range(100., 1500., as_db=False)
print("Selected Database has %g records" % len(db1))
# Eliminate records missing Rrup
selector2 = SMRecordSelector(db1)
db1.records = selector2.select_within_distance_range("rrup", 0.1, 500., alternative=False,as_db=False)
print("Selected Database has %g records" % len(db1))
Now select the records only within a bounding box around Italy (defined by west, south, east, north)
In [ ]:
# Italy Bounding Box
west, south, east, north = (5.0, 35.0, 20., 47.)
# Select the records
selector3 = SMRecordSelector(db1)
db1.records = selector2.select_epicentre_within_bounding_box(west, south, east, north, as_db=False)
print("Selected Database has %g records" % len(db1))
# Look at the magnitude-distance distribution
db_magnitude_distance(db1, "rrup")
In [ ]:
resid1 = res.Residuals(gsim_list, imt_list)
resid1.get_residuals(db1, component="SARotD50")
print("Printing to file")
resid1.pretty_print("NGAWest2Residuals1.csv", sep=",")
Explore the overall residual trend
In [ ]:
rspl.ResidualPlot(resid1, "BooreEtAl2014", "PGA", figure_size=(8,8))
rspl.ResidualPlot(resid1, "BindiEtAl2014Rjb", "PGA", figure_size=(8,8))
rspl.ResidualPlot(resid1, "CampbellBozorgnia2014", "PGA", figure_size=(8,8))
rspl.ResidualPlot(resid1, "CauzziEtAl2014", "PGA", figure_size=(8,8))
In [ ]:
rspl.ResidualPlot(resid1, "BooreEtAl2014", "SA(1.0)", figure_size=(8,8))
rspl.ResidualPlot(resid1, "BindiEtAl2014Rjb", "SA(1.0)", figure_size=(8,8))
rspl.ResidualPlot(resid1, "CampbellBozorgnia2014", "SA(1.0)", figure_size=(8,8))
rspl.ResidualPlot(resid1, "CauzziEtAl2014", "SA(1.0)", figure_size=(8,8))
Compare with Magnitude
In [ ]:
rspl.ResidualWithMagnitude(resid1, "BooreEtAl2014", "PGA", figure_size=(8,8))
rspl.ResidualWithMagnitude(resid1, "BindiEtAl2014Rjb", "PGA", figure_size=(8,8))
rspl.ResidualWithMagnitude(resid1, "CampbellBozorgnia2014", "PGA", figure_size=(8,8))
rspl.ResidualWithMagnitude(resid1, "CauzziEtAl2014", "PGA", figure_size=(8,8))
In [ ]:
rspl.ResidualWithMagnitude(resid1, "BooreEtAl2014", "SA(1.0)", figure_size=(8,8))
rspl.ResidualWithMagnitude(resid1, "BindiEtAl2014Rjb", "SA(1.0)", figure_size=(8,8))
rspl.ResidualWithMagnitude(resid1, "CampbellBozorgnia2014", "SA(1.0)", figure_size=(8,8))
rspl.ResidualWithMagnitude(resid1, "CauzziEtAl2014", "SA(1.0)", figure_size=(8,8))
Compare with Distance
In [ ]:
rspl.ResidualWithDistance(resid1, "BooreEtAl2014", "SA(1.0)", figure_size=(8,8), plot_type="linear")
rspl.ResidualWithDistance(resid1, "BindiEtAl2014Rjb", "SA(1.0)", figure_size=(8,8), plot_type="linear")
rspl.ResidualWithDistance(resid1, "CampbellBozorgnia2014", "SA(1.0)", figure_size=(8,8), plot_type="linear")
rspl.ResidualWithDistance(resid1, "CauzziEtAl2014", "SA(1.0)", figure_size=(8,8), plot_type="linear")
Run the Log-likelihood Analysis and see Results: 1. By Period, 2. By GMPE
In [ ]:
llh = res.LLH(gsim_list, imt_list)
llh.get_residuals(db1, component="SARotD50")
In [ ]:
# Run LLH analysis and store the results in dictionaries
llh_values, weights = llh.get_loglikelihood_values(imt_list)
periods = dict([(gmpe, []) for gmpe in gsim_list])
llhs = dict([(gmpe, []) for gmpe in gsim_list])
for gmpe in gsim_list:
print("LLH(%s) = %12.6f Weight = %12.6f" % (gmpe, llh_values[gmpe]["All"], weights[gmpe]))
for imt in imt_list:
if llh_values[gmpe][imt]:
print(" %s = %12.6f" % (imt, llh_values[gmpe][imt]))
llhs[gmpe].append(llh_values[gmpe][imt])
if "SA(" in imt:
periods[gmpe].append(from_string(imt).period)
elif "PGA" in imt:
periods[gmpe].append(0.0)
else:
pass
# Plot the LLH values by Period
fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)
for gmpe in gsim_list:
ax.plot(periods[gmpe], llhs[gmpe], marker="o", lw=1.5,
label="{:s} ({:.4f})".format(gmpe, llh_values[gmpe]["All"]))
ax.set_xlabel("Period (s)", fontsize=14)
ax.set_ylabel("LLH", fontsize=14)
ax.grid(True)
ax.legend(fontsize=14)
In [ ]: